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Abstract 

The inositol trisphosphate receptor (IP3R) is one of the most important cellular components responsible for oscillations in 
the cytoplasmic calcium concentration. Over the past decade, two major questions about the IP3R have arisen. Firstly, how 
best should the IP 3 R be modeled? In other words, what fundamental properties of the IP 3 R allow it to perform its function, 
and what are their quantitative properties? Secondly, although calcium oscillations are caused by the stochastic opening 
and closing of small numbers of IP3R, is it possible for a deterministic model to be a reliable predictor of calcium behavior? 
Here, we answer these two questions, using airway smooth muscle cells (ASMC) as a specific example. Firstly, we show that 
periodic calcium waves in ASMC, as well as the statistics of calcium puffs in other cell types, can be quantitatively 
reproduced by a two-state model of the IP3R, and thus the behavior of the IP 3 R is essentially determined by its modal 
structure. The structure within each mode is irrelevant for function. Secondly, we show that, although calcium waves in 
ASMC are generated by a stochastic mechanism, IP3R stochasticity is not essential for a qualitative prediction of how 
oscillation frequency depends on model parameters, and thus deterministic IP 3 R models demonstrate the same level of 
predictive capability as do stochastic models. We conclude that, firstly, calcium dynamics can be accurately modeled using 
simplified IP3R models, and, secondly, to obtain qualitative predictions of how oscillation frequency depends on 
parameters it is sufficient to use a deterministic model. 
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Introduction 

Oscillations in cytoplasmic calcium concentration ([Ca 2+ ]j), 
mediated by inositol trisphosphate receptors (IP3R; a calcium 
channel that releases calcium ions (Ca 2 + ) from the endoplasmic or 
sarcoplasmic reticulum (ER or SR) in the presence of inositol 
trisphosphate (IP3)) play an important role in cellular function in 
many cell types. Hence, a thorough knowledge of the behavior of 
the IP3R is a necessary prerequisite for an understanding of 
intracellular Ca 2+ oscillations and waves. Mathematical and 
computational models of the IP3R play a vital role in studies of 
Ca 2+ dynamics. However, over the past decade, two major 
questions about IP3R models have arisen. 

Firstly, how best should the IP3R be modeled? Models of the 
IP3R have a long history, beginning with the heuristic models of [1— 
3] . With the recent appearance of single-channel data from IP3R in 
vivo [4,5], a new generation of Markov IP 3 R models has recently 
appeared [6,7]. These models show that IP3R exist in different 
modes with different open probabilities. Within each mode there 
are multiple states, some open, some closed. Importandy, it was 
found [8] that time-dependent transitions between different modes 



are crucial for reproducing Ca 2 + puff data from [9] . However, it is 
not yet clear whether transitions between states within each mode 
are important, or whether all the important behaviors are captured 
simply by inter-mode transitions. 

Secondly, why do deterministic models of the IP3R perform so 
well as predictive models? Deterministic models of the IP3R have 
proven to be useful predictive models in a range of cell types. For 
example, IPsR-based models have been developed to study Ca 2+ 
oscillations in airway smooth muscle cells (ASMC) [10-13], and 
these models have made predictions which have been confirmed 
experimentally. This shows the usefulness of such models in 
advancing our understanding of how intracellular Ca 2+ oscilla- 
tions and waves are initiated and controlled in ASMC. However, 
these models are deterministic models which assume infinitely 
many IP3R per unit cell volume, an assumption that contradicts 
experimental findings in many cell types showing that Ca 2+ puffs 
and spikes occur stochastically, and that intracellular Ca 2+ waves 
and oscillations arise as an emergent property of fundamental 
stochastic events [9,14,15]. 

Here, we answer these two fundamental modeling questions 
using data and models from ASMC. Firstly, we show that a simple 
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Author Summary 

The inositol trisphosphate receptor (IP3R) is one of the 
most important cellular components responsible for 
calcium oscillations. Over the past decade, two major 
questions about the IP 3 R have arisen. Firstly, what 
fundamental properties of the IP 3 R allow it to perform 
its function? Secondly, although calcium oscillations are 
caused by the stochastic properties of small numbers of 
IP3R is it possible for a deterministic model to be a reliable 
predictor of calcium dynamics? Using airway smooth 
muscle cells as an example, we show that calcium 
dynamics can be accurately modeled using simplified 
IP3R models, and, secondly, that deterministic models are 
qualitatively accurate predictors of calcium dynamics. 
These results are important for the study of calcium 
dynamics in many cell types. 



model of the IP3R, involving only two states with time-dependent 
transitions, suffices to generate correct dynamics of Ca 2+ puffs 
and oscillations. Secondly, we show that, although Ca 2+ 
oscillations in ASMC are generated by a stochastic mechanism, 
a deterministic model can make the same qualitative predictions as 
the analogous stochastic model, indicating that deterministic 
models, that require much less computational time and complex- 
ity, can be used to make reliable predictions. Although we work in 
the specific context of ASMC, our results are applicable to other 
cell types that exhibit similar Ca 2+ oscillations and waves. 

Results 

A two-state model of the IP3R is sufficient to reproduce 
function 

We have previously shown [8] that the statistics of Ca 2+ puffs in 
SH-SY5Y cells can be reproduced by a Markov model of the IP3R 
based on the steady-state data of [5] and the time-dependent data 
of [4]. In this model the IP3R can exist in 6 different states, 
grouped into two modes, which we call Drive and Park (see 
Fig. 1). The Drive mode (which contains 4 states; 1 open and 3 
closed) has an average open probability of around 0.7, while the 
Park mode (which contains the remaining two states; 1 open and 1 
closed) has an open probability close to zero. Transitions between 
states within each mode are independent of Ca 2+ and IP3; only 
the transitions between modes are ligand-dependent. 

In our previous study on calcium puffs [8], we showed that, to 
reproduce the experimentally observed non-exponential interspike 
interval (ISI) distribution and coefficient of variation (CV) of ISI 
smaller than 1, the time-dependent intermodal transitions are 
crucial. Lack of time dependencies in the Siekmann model leads to 
exponential ISI distributions and C V = 1 , which is not the case for 
calcium spikes in ASMC. Fig. 2 A shows an example of Ca 2+ 
oscillations generated by 50 nM methacholine (MCh, an agonist 
that can induce the production of IP3 by binding to a G protein- 
coupled receptor in the cell membrane) in ASMC. By gathering 
data from 14 cells in 5 mouse lung slices, we found that the 
standard deviation of the interspike interval (ISI) is approximately 
a linear function of the ISI mean, with a slope clearly between 0 
and 1 (i.e. CV< 1), indicating that the spikes are generated by an 
inhomogeneous Poisson process (a slope of 1 would denote a pure 
Poisson process) (see Fig. 2B). This shows the necessity of inclusion 
of time-dependent transitions for mode-switching. 

Using a quasi-steady-state approximation, and ignoring states 
with very low dwell times, it is possible to construct a simplified 
two-state version of the full six-state model (see Materials and 




Figure 1. The structure of the Siekmann IP3R model. The IP 3 R 
model is comprised of two modes. One is the drive mode containing 
three closed states C\, C%, C3 and one open state 0$. The other is the 
park mode which includes one closed state C4 and one open state O5. 
qs are rates of state-transitions between two adjacent states and 942 
and §24 are transitions between the two modes [7]. 
doi:10.1371/journal.pcbi.1003783.g001 

Methods). In the simplified model the intramodal structure is 
ignored, and only the intermodal transitions have an effect on 
IP3R behavior. In Fig. 3 we compared the simplified IP3R model 
to the full six-state model. Both models have the same distribution 
of interspike interval, spike amplitude and spike duration. 
Moreover, by looking at a more detailed comparison between 
the two model results (Figs. 4A, C and E) and experimental data 
(Figs. 4B, D and F), we found the 2-state model not only can 
reproduce the behaviour of the 6-state model, but can also 
qualitatively reproduce experimental data. The average experi- 
mental ISI shows a clear decreasing trend as MCh concentration 
increases (although a saturation occurs in the data for high MCh), 
a trend that is mirrored by the model results as the IP3 
concentration increases. Unfortunately, since the exact relation- 
ship between MCh concentration and IP3 concentration is 
uncertain, a quantitative comparison is not possible. In both 
model and experimental results, the average peak and duration of 
the oscillations are nearly independent of agonist concentration. 
The quantitative difference in spike duration between the model 
results and the data in Figs. 4E and F are most likely due to choice 
of calcium buffering parameters. For example, adding 3 or 5 fiM 
fast Ca 2+ buffer (see Materials and Methods) increases the average 
spike duration to 0.54 s or 0.7 s respectively, which are close to the 
levels shown in the data. 

Thus, the intramodal structure of the six-state model is 
essentially unimportant, as the model behavior (in terms of the 
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time(s) mean of ISIs (s) 

Figure 2. Ca 2 oscillations in ASMC in lung slices are generated by a stochastic mechanism. A: experimental Ca 2+ spiking in ASMC in 
lung slices, stimulated with 50 nM MCh. In the upper panel we filter out baseline noise by using a low threshold of 1.42 (relative fluorescence 
intensity) and then choose samples with amplitude larger than 1.75. The ISI calculated from the upper panel is shown in the lower panel. B: 
relationship between the standard deviation and the mean of experimental ISIs. Data obtained from 1 4 ASMC in 5 mouse lung slices. The relationship 
is approximately linear with a slope of 0.66, which implies that an inhomogeneous Poisson process governs the generation of oscillations. The dashed 
line indicates where the coefficient of variation (CV) is 1 (as it is for a pure Poisson process). Variation in ISI is mainly caused by both use of different 
doses of MCh and different sensitivities of different cells to MCh. Error bars indicate the standard errors of the means (SEM). 
doi:10.1371/journal.pcbi.1003783.g002 



statistics of puffs and oscillations) is governed almost entirely by the 
time dependence of the intermode transitions, particularly the time 
dependence of the rapid inhibition of the IP3R by high [Ca ];, 
and the slow recovery from inhibition by Ca 2+ . The multiple 
states within each mode are necessary to obtain an acceptable 
quantitative fit to single-channel data, but are nevertheless of 
limited importance for function. Hence, even when simulating 
microscopic events such as Ca 2+ puffs it is sufficient to use a 
simpler, faster, two-state model, rather than a more complex six- 



state model. In the following, we will use the 2-state IP3R model to 
generate all the simulation results. 

Prediction of stochastic Ca 2+ behavior by a deterministic 
model 

Although the data (Fig. 2) show that Ca 2+ oscillations in ASMC 
are generated by a stochastic process, not a deterministic one, we 
wish to know to what extent a deterministic model can be used to 
make qualitative (and experimentally testable) predictions. Our 
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Figure 3. A 2-state open/closed model quantitatively reproduces the 6-state IP3R model. A: histograms of interspike interval (ISI) 
distribution for both the 6-state and the simplified models. The ISI is defined to be the waiting time between successive spikes. Each histogram 
contain an equal number of samples (180). B: comparison of average ISI, average peak value of [Ca 2+ ]; (c in the model) and average spike duration. 
All distributions were computed at a constant [IP 3 ] = 0.15 jiM. 
doi:10.1371/journal.pcbi.1003783.g003 
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Figure 4. More detailed comparisons between the 2-state and the 6-state IP 3 R models, and a comparison to experimental data. As a 

function of IP 3 concentration (p), the two models give the same ISI (A), peak [Ca 2+ ]j (C) and spike duration (E). These results agree qualitatively with 
experimental data, as shown in panels B, D and F respectively. Quantitative comparisons are generally not possible as the relationship between IP 3 
concentration and agonist concentration is not known. Error bars represent mean + SEM. Data for each MCh concentration are obtained from at 
least three different cells from at least two different lung slices. 
doi:10.1371/journal.pcbi.1003783.g004 
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simplified 2-state Markov model of the IP3R can be converted to a 
deterministic model (see Materials and Methods). The result is a 
system of ordinary differential equations (ODEs) with four 
variables, which takes into account the increased [Ca 2+ ] ; at an 
open IP3R pore, as well as the increased [Ca 2 + ]j within a cluster 
of IP3R; the four variables are the [Ca 2+ ] ; outside the IP 3 R 
cluster (c), the [Ca 2+ ]; within the IP 3 R cluster (q,), the total 
intracellular Ca 2+ concentration (c t ) and an IP3R gating variable 
(A42). We refer to the reduced 4D model as the deterministic model 
for all the results and analyses. 

Note that there is no physical or geometric constraint enforcing 
a high local [Ca 2+ ] ; ; in this case the spatial heterogeneity arises 
solely from the low diffusion coefficient of Ca 2+ . Our use of cj is 
merely a highly simplified way of introducing spatial heterogeneity 
of the Ca 2+ concentration. Since the IP3R can only "see" c/, (as 
well as the Ca 2+ concentration right at the mouth of an open 
channel, which we denote by c p ), but cannot be influenced directly 
by c (the experimentally observed Ca 2+ signal), our approach 
allows for the functional differentiation of the rapid local 
oscillatory Ca 2+ in the cluster, from the slower Ca 2+ signal in 
the cytoplasm, without the need for computationally intensive 
simulations of a partial differential equation model. Quantitative 
accuracy is thus sacrificed for computational convenience. 

Calcium oscillations in the stochastic and deterministic models 
are shown in Fig. 5A. According to our previous results [8], the 
average value of A42 over the cluster of IP3R primarily regulates 
the termination and regeneration of individual spikes. This can be 
seen in the stochastic model by projecting the solution on the 
c b, ^42 phase plane (Fig. 5B). Upon an initial Ca 2+ release from 
one or more IP3R, a large spike is generated by Ca 2+ -induced 
Ca 2+ release (via the IP3R) during which time a decreasing A42 
gradually decreases the average open probability of the clustered 
IP3R. The spike is terminated when A42 is too small to allow 
further Ca 2+ release. This phenomenon is qualitatively repro- 
duced by the deterministic model (Fig. 5D). In both the stochastic 
and deterministic models the decrease in average IP3R open 
probability of a cluster of IP3R caused by Ca 2+ inhibition is the 
main reason for the termination of each spike. 

According to Figs. 5B and D, regeneration of each spike 
requires a return of ^42 back to a relatively high value (i.e., 
recovery of the IP3R from inhibition by Ca 2+ ). The deterministic 
model sets a clear threshold for the regeneration, as can be seen in 
Fig. 5C, where an upstroke in cj occurs when the trajectory creeps 
beyond the sharp "knee" of the white curve. When the trajectory 
reaches the knees of the white curve it is forced to jump across to 
the other stable branch of the critical manifold, resulting in a fast 
increase in c/, followed by a relatively fast increase in c (seen by 
combining Figs. 5C and D). 

In contrast, the stochastic model enlarges the contributions of 
individual IP3R so that the generation of each spike is also 
effectively driven by random Ca 2+ release through the IP3R, 
which can be seen in the inset of Fig. 5B where the site of spike 
initiation (blue bar) exhibits significantly greater variation than 
that of spike termination (green bar). In spite of this, the essential 
similarities in phase plane behavior result in both deterministic 
and stochastic models making the same qualitative predictions in 
response to perturbations, such as changes in IP3 concentration 
([IP3]), Ca 2+ influx or efflux. In the following, we illustrate this by 
investigating a number of experimentally testable predictions. Due 
to the extensive importance of frequency encoding in many Ca 2+ - 
dependent processes, we focus particularly on the change of 
oscillation frequency in response to parameter perturbations. As a 



side issue we also investigate how the oscillation baseline depends 
on physiologically important parameters. 

Dependence of oscillation frequency on IP 3 concentration 

In many cell types a moderate increase in [IP3] increases the 
Ca 2+ oscillation frequency (see Fig. 2A in [11], Fig. 4E in [16] 
and Fig. 6B in [17]), a result that is reproduced by both model 
types (Fig. 6A). As [IP3] increases, the stochastic model increases 
the probability of the initial Ca 2+ release through the first open 
IP3R and of the following Ca 2+ release, thus shortening the 
average ISI. Although the oscillatory region of the deterministic 
model is strictly confined by bifurcations which do not apply to the 
stochastic model, the deterministic model can successfully replicate 
an increasing frequency by lowering the "knee" of the red curve in 
Fig. 5D and shortening the time spent from the termination point 
c to the initiation point a (thus shortening the ISI). Hence, 
although the deterministic model cannot be used to predict the 
exact values of [IP3] at which the oscillations begin and end, as 
stochastic effects predominate in these regions, it can be used to 
predict the correct qualitative trend in oscillation frequency. 

Dependence of oscillation frequency on Ca 2+ influx and 
efflux 

In many cell types, including ASMC, transmembrane fluxes 
modulate the total intracellular Ca 2+ load (c t ) on a slow time scale 
[16,18], and thereby modulate the oscillation frequency [19]. 
Experimental data can be seen in Fig. 8 in [16] and Fig. 2 in [18]. 
Figs. 6B and C show that both stochastic and deterministic models 
predict the same qualitative changes in oscillation frequency in 
response to changes in membrane fluxes (through membrane 
ATPase pumps and/ or Ca 2+ influx channels such as receptor- 
operated channels or store-operated channels). 

Dependence of oscillation frequency on SERCA 
expression 

The level of sarco/endoplasmic reticulum calcium ATPase 
(SERCA) expression (or capacity) is important for airway 
remodeling in asthma [20] and ASMC Ca 2+ oscillations [21]. 
We thus investigated the predictions of the two models in response 
to changes in SERCA expression (V s ). As V s decreases, the 
deterministic model exhibits a decreasing frequency, in agreement 
with experimental data (see Figs. 3 and 4 in [21]). The same trend 
is seen in the stochastic model with only 20 IP3R (see Fig. 6D). 

Dependence of oscillation frequency on Ca 2+ buffer 
concentration 

Calcium buffers have been shown to be able to change the ISI 
and spike duration, which in turn change the oscillation frequency 
[15,22]. We compared the effects on the two models of varying 
total buffer concentration (B t ) by adding one buffer with relatively 
fast kinetics to the models (see Materials and Methods for details). 
In both models the frequency decreases as B t increases (see 
Fig. 6E), which is consistent with experimental data (Fig. 2B in 
[18]). This is not surprising, because increasing B t can decrease 
the effective rates of SR Ca 2+ release and reuptake. 

Dependence of oscillation baseline on Ca 2+ influx and 
SERCA expression 

Sustained elevations of baseline during agonist- induced Ca 2+ 
oscillations or transients have been observed experimentally, and 
are believed to be a result of an increase in Ca 2+ influx caused by 
opening of membrane Ca 2+ channels [13,16]. Furthermore, there 
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Figure 5. Stochastic and deterministic simulations exhibit similar dynamic properties. A: simulated stochastic (upper panel) or 
deterministic (lower panel) Ca 2+ oscillations at 0.1 /iM IP3. B: a typical stochastic solution projected on the ej— A42 plane. The average A42 
represents the average value of A42 over the 20 IP3R. Statistics (mean + SD) of the initiation point (blue square), the peak (red square) and 
termination point (green square) are shown in the inset. 1 16 samples are obtained by applying a low threshold of 0.15 [iM and a high threshold of 
0.8 /iM to a. C: a typical periodic solution of the deterministic model (black curve), plotted in the c, cj, I142 phase space. The arrow indicates the 
direction of movement, c, is the slowest variable so that its variation during an oscillation is very small. This allows to treat c, as a constant 
(c, = 53.12 [iM in this case) and study the dynamics of the model in the c, c*, I142 phase space. The color surface is the surface where tfej/rfi=0 (called 
the critical manifold). The white N-shaped curve is the intersection of the critical manifold and the surface dc/dt = 0. D: projection of the periodic 
solution to the cj, ha plane. The red N-shaped curve is the projection to the cj, plane of the white curve shown in C. The evolution of the 
deterministic solution exhibits three different time scales separated by green circles (labelled by a, b and c) and indicated by arrows (triple arrow: 
fastest; double arrow: intermediate; single arrow: slowest). 
doi:10.1371/journal.pcbi.1003783.g005 



is evidence showing that decreased SERCA expression could also 
increase the baseline (Fig. 4 in [21]). Those phenomena are 
successfully reproduced by both models (see Fig. 7). 

Discussion 

In this paper we address two current major questions in the field 
of Ca 2+ modeling. Firstly, we show that Ca 2+ puffs and stochastic 
oscillations can be reproduced quantitatively by an extremely 
simple model, consisting only of two states (one open, one closed), 
with time-dependent transitions between them. This model is 
obtained by removing the intramodal structure of a more complex 
model that was determined by fitting a Markov model to single- 



channel data [7] . We thus show that the internal structure of each 
mode is irrelevant for function and mode switching is the key 
mechanism for the control of calcium release. The necessity for 
time-dependent mode switching is shown not only by the dynamic 
single-channel data of [4]), but also by the puff data of [9] and our 
ASMC data. 

Secondly, we investigate the role of stochasticity of IP3R in 
modeling Ca 2 + oscillations in ASMC by comparing a stochastic 
IP 3 R-based Ca 2+ model and its associated deterministic version, 
for parameters such that both of the models exhibit Ca 2 + spikes 
but the stochastic model cannot necessarily be replaced by a 
mean-field model. We find that a four-variable deterministic 
model has the same predictive power as the stochastic model, in 
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Figure 6. Comparison of parameter-dependent frequency changes in the stochastic and deterministic models. All curves are 
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doi:10.1371/journal.pcbi.1003783.g006 
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Figure 7. Dependence of calcium oscillation baseline on calcium influx and SERCA expression. A: increasing influx (described by V mn ) 
increases the average trough of Ca 2+ oscillations. B: decreasing SERCA expression (described by V s ) increases the average trough of Ca 2+ 
oscillations. All curves are computed at 0.12 /iM IP3. 
doi:10.1371/journal.pcbi.1003783.g007 



that it correctly reproduces the process of spike termination and 
predicts the same qualitative changes in oscillation frequency and 
baseline in response to a variety of perturbations that are 
commonly used experimentally. The mechanism for termination 
of individual spikes is fundamentally a deterministic process 
controlled by a rapid inhibition induced by the high local [Ca 2+ ] ; 
in the IP3R cluster, whereas spike initiation is significantly affected 
by stochastic opening of IP3R. Hence, repetitive Ca 2+ cycling is 
primarily induced by the time-dependent gating variables govern- 
ing transitions of the IP3R from one mode to another. 

Our simplified two-state model of the IP3R is identical in structure 
(although not in parameter values) to the well-known model of [23] . It 
is somewhat ironic that after 20 years of detailed studies of the IP3R 
and the construction of a plethora of models of varying complexity, 
the single-channel data have led us around full circle, back to these 
original formulations. Excitability is arising via a fast activation 
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Figure 8. Schematic diagram of the Ca 2+ model, c represents 
cytoplasmic Ca 2+ concentration, excluding a small local Ca 2+ (whose 
concentration is denoted by q,) close to the Ca 2+ release site (i.e., an 
IP3R cluster). Upon coordinated openings of the IP3R, SR Ca 2+ (cj is 
first released into the local domain (Jipr) to cause a rapid increase in q,. 
High local Ca 2+ then diffuses to the rest of the cytoplasm (/dtff)/ and is 
eventually pumped back to the SR (/ scrca ). 
doi:10.1371/journal.pcbi.1003783.g008 



followed by a slower inactivation, a combination often seen in 
physiological processes [24] . Encoding of this fundamental combi- 
nation results directly from the two-mode structure of the IP3R. 
Although similar single-channel data have been used to construct 
three-mode models [6,25], neither of these models has yet been used 
in detailed studies of Ca 2+ puffs and waves, and it remains unclear 
whether or not they have a similar underlying structure. 

In contrast to previous deterministic ODE models, our four- 
variable Ca 2+ model includes a more accurate IP3R model, as 
well as local control of clustered IP3R by two distinct Ca 2+ 
microdomains; one at the mouth of an open IP3R, the other inside 
a cluster of IP3R. Neglect of either of these microdomains leads to 
models that either exhibit unphysiological cytoplasmic Ca 2+ 
concentrations or fail to reproduce reasonable oscillations. This 
underlines the importance of taking Ca 2+ microdomains into 
consideration when constructing any model. Our microdomain 
model is highly simplified, with the microdomain being treated 
simply as a well-mixed compartment. More detailed modeling of 
spatially-dependent microdomains is possible, and not difficult in 
principle, but requires far greater computational resources. It is 
undeniable that a more detailed model, incorporating the full 
spatial complexity - and possibly stochastic aspects as well - would 
make, overall, a better predictive tool. However, our goal is to find 
the simplest models that can be used as predictive tools. 

An important similar study is that of Shuai and Jung [26] . They 
compared the use of Markov and Langevin approaches to the 
computation of puff amplitude distributions, compared their results 
with file deterministic limit, and showed that IP3R stochasticity does 
not qualitatively change the type of puff amplitude distribution except 
for when there are fewer than 10 IP3R. Here, we significantly extend 
the scope of their study by exploring the effects of IP3R stochasticity on 
the dynamics of Ca 2+ spikes, and we do this in the context of an IP3R 
model that has been fitted to single-channel data. Although this is true 
in a general sense for the Li-Rinzel model, which is based on the 
DeYoung-Keizer model, which did take into account the opening time 
distributions of IP3R in lipid bilayers, neither model can reproduce the 
more recent data obtained from on-nuclei patch clamping. When these 
recent data are taken into account one obtains a model with the same 
structure, but quite different parameters and behavior. 
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Table 1. Parameter values of the stochastic calcium model. 



Parameter 


Description 


Value/Units 




IP3R flux coefficient 


0.05 s" 1 


fcdiff 


Ca 2+ diffusional flux coefficient 


10s-' 


'''leak 


SR leak flux coefficient 


0.0032 s" 1 


V, 


maximum capacity of SERCA 


10 /jM-s- 1 


K„ 


SERCA half-maximal activating [Ca 2+ ]; 


0.26 


n s 


Hill coefficient for SERCA 


1.75 




plasma membrane leak influx 


0.03115 fiMs-' 


K rocc 


ROCC flux coefficient 


0.2 s~' 


^SOCC 


maximum capacity of SOCC 


1.6 fiMs- 1 




SOCC dissociation constant 


100 t M 


y, 


maximum capacity of plasma pump 


0.8 ((M-s- 1 




half-maximal activating [Ca 2 + ] i of plasma pump 


0.5 fiM 


"p 


Hill coefficient for plasma pump 


2 


Yl 


the cytoplasmic-to-microdomain volume ratio 


100 


h 


the cytoplasmic-to-SR volume ratio 


10 




an instantaneous high [Ca 2 + ] i at open channel 
pore when c s = 100 


120 t M 


N, 


total number of IP 3 R channels 


20 



doi:l 0.1 371 /journal.pcbi.l 003783.t001 



We find that, in spite of a relatively large variation in spike 
amplitude which is partially caused by a large variation in ISI 
(Fig. 5B), the mechanism governing individual spike terminations 
is the same for both a few or infinitely many IP3R, which explains 
why the one-peak type of amplitude distribution is independent of 
the choice of IP3R number (see Fig. 6A in [26]). 

Another important relevant study was done by Dupont et al. [27] , 
who compared the regularity of stochastic oscillations in hepatocytes 
for different numbers of IP3R clusters. They found that the impact 
of IP3R stochasticity on global Ca 2 + oscillations (in terms of CV) 
increases as the total cluster number decreases. Our study here 
extends these results, and demonstrates how well stochastic 
oscillations can be qualitatively described by a deterministic system, 
even when there is only a small number of IP3R (which appears to 
be the case for ASMC, in which the wave initiation site is only 
2 ~ 4 /mi in diameter). Indeed, as we have shown, for the purposes 
of predictive modeling a simple deterministic model does as well as 
more complex stochastic simulations. 

Ryanodine receptors (RyR) are another important component 
modulating ASMC Ca 2+ oscillations [16,28,29] but are not 
included in our model. This is because the role of RyR is not fully 
understood and may be species-dependent; for example, in mouse 
or human ASMC, RyR play very little role in IP3 -induced 
continuing Ca 2+ oscillations [17,30], but this appears not to be 
true for pigs [28]. Our study focuses on the calcium oscillations in 
mouse and human (as we did in our experiments) where inclusion 
of a deterministic model of RyR should have little effect. An 
understanding of the role of RyR stochasticity and how the IP3R 
and the RyR interact needs a reliable RyR Markov model, 
exclusive to ASMC, which is not currently available. Multiple 
Markov models of the RyR have been developed for use in cardiac 
cells [31], but these are based on single-channel data from lipid 
bilayers, and are adapted for the specific context of cardiac cells. 
Their applicability to ASMC remains unclear. 



Although we have not shown that the deterministic model 
for ASMC has the same predictive power as the stochastic 
model in all possible cases (which would hardly be possible in 
the absence of an analytical proof) the underlying similarity in 
phase plane structure indicates that such similarity is plausible at 
least. Certainly, we have not found any counterexample to this 
claim. However, whether or not this claim is true for all cell 
types is unclear. Some cell types exhibit both local Ca 2+ puffs 
and global Ca 2+ spikes (usually propagating throughout the 
cells in the form of traveling waves), showing that initiation of 
such Ca 2+ spikes requires a synchronization of Ca 2+ release 
from more than one cluster of IP3R [14]. This type of spiking 
relies on the hierarchical organization of Ca 2+ signal pathways, 
in particular the stochastic recruitment of both individual IP3R 
and puffs at different levels [32], and therefore cannot be simply 
reproduced by deterministic models containing only a few 
ODEs. However, Ca 2+ oscillations in ASMC, as observed in 
lung slices, may not be of this type, as IP 3 R-dependent puffs 
have not been seen in these ASMC. It thus appears that, in 
ASMC in lung slices, every Ca 2+ "puff" initiates a wave, 
resulting in periodic waves with ISI that are governed by the 
dynamics of individual puffs. 

Materials and Methods 

Ethics Statement 

Animal experimentations carried out were approved by the 
Animal Care and Use Committee of the University of Massachu- 
setts Medical School under approval number A-836-12. 

Lung slice preparation 

BALB/c mice (7-10 weeks old, Charles River Breeding Labs, 
Needham, MA) were euthanized via intraperitoneal injection of 
0.3 ml sodium pentabarbitone (Oak Pharmaceuticals, Lake Forest, 
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IL). After removal of the chest wall, lungs were inflated with 
~1.1 ml of 1.8% warm agarose in sHBSS via an intratracheal 
catheter. Subsequently, air (~ 0.3 ml) was injected to push the 
agarose within the airways into the alveoli. The agarose was 
polymerized by cooling to 4 a C. A vibratome (VF-300, Precisionary 
Instruments, San Jose, CA) was used to make 180 /im thick slices 
which were maintained in Dulbecco's Modified Eagle's Media 
(DMEM, Invitrogen, Carlsbad, CA) at 37° C in 10% C0 2 /air. All 
experiments were conducted at 37°C in a custom-made temper- 
ature-controlled Plexiglas chamber as described in [17]. 

Measurement of Ca 2+ oscillations 

Lung slices were incubated in sHBSS containing 20 pM 
Oregon Green 488 BAPTA-1-AM (Invitrogen), a Ca 2+ -indicator 
dye, 0.1% Pluronic F-127 (Invitrogen) and 200 /iM sulfobromo- 
phthalein (Sigma Aldrich, St Louis, MO) in the dark at 30°C for 
1 hour. Subsequently, the slices were incubated in 200 fiM 
sulfobromophthalein for 30 minutes. Slices were mounted on a 
cover-glass and held down with 200 fim mesh. A smaller cover- 
glass was placed on top of the mesh and sealed at the sides with 
silicone grease to facilitate solution exchange. Slices were 
examined with a custom-built 2-photon scanning laser microscope 
with a x 40 oil immersion objective lens and images recorded at 
30 images per second using Videosavant 4.0 software (IO 
Industries, Montreal, Canada). Changes in fluorescence intensity 
(which represent changes in [Ca 2+ ];) were analyzed in an ASMC 
of interest by averaging the grey value of a 1 0 x 10 pixel region 
using custom written software. Relative fluorescence intensity 
(F/Fo) was expressed as a ratio of the fluorescence intensity at a 
particular time (F) normalized to the initial fluorescence intensity 
(Fo). 

The calcium model 

Inhomogeneity of cytoplasmic Ca 2+ concentration not only 
exists around individual channel pores of the IP3R, where a nearly 
instantaneous high Ca 2+ concentration at the pore (denoted by c p ) 
leads to a very sharp concentration profile, but is also seen inside 
an IP3R cluster where the average cluster Ca 2+ concentration (ct) 
is apparently higher than that of the surrounding cytoplasm (c) 
[33]. This indicates that during Ca 2+ oscillations each IP3R is 
controlled by either the pore Ca 2 + concentration (when it is open) 
or the cluster Ca 2+ concentration (when it is closed). Neither of 
these local concentrations influence cell membrane fluxes or the 
majority of SERCAs, which we assume to be distributed outside 
the cluster. 

The scale separation between the pore Ca 2+ concentration and 
the cluster Ca 2+ concentration allows to treat c p as a parameter, 
providing a simpler way of modeling local Ca 2+ events (like Ca 2+ 
puffs) that has been used in several previous studies [8,34,35]. 
However, evolution of the cluster concentration and wide-field 
cytoplasm Ca 2+ concentration are not always separable, so an 
additional differential equation for the cluster Ca 2+ is necessary. 

A schematic diagram of the model is shown in Fig. 8. The 
corresponding ODEs are 

dc 

~~T = /dif f 4~ /leak J serca 4" Jin J pm > ( 1 ) 



~jT =yi(^n>R— ^diff)> (2) 



i. — /in /prm (3) 

where c t = c + C[,/y i -\-c s /y 2 representing total intracellular Ca 2+ 
concentration, and thus SR Ca 2+ concentration, c s is given by 
c s = y 2 (c, — c — ci,/y- l ). y^ and y 2 are the volume ratios given in 
Table 1. /ipr is the flux through the IP3R, /i ea k is a background 
Ca 2+ leak out of the SR, and / se ica is the uptake of Ca 2+ into the 
SR by SERCA pumps. / pm is the flux through plasma pump, and 
/in represents a sum of main Ca 2+ influxes including / roC c 
(receptor-operated Ca 2+ channel), / SOC c (store-operated Ca 2+ 
channel) and /leakin (Ca 2+ leak into the cell). Jdiff coarsely models 
the diffusion flux from cluster microdomain to the cytoplasm. 
Details of the fluxes are 

• Different formulations of /ipr give different types of models: 

a) For the stochastic model, /ipr = (kipjt./Nt)N 0 (c s — c) where 
^ipr is the maximum conductance of a cluster of N t IP3R 
(here N t = 20). N 0 is the number of open IP3R determined by 
the states of IP 3 R. 

b) For the deterministic model we set /ipr = /:tpr P 0 (c s — c) 
where P„ is the IP3R open probability, a continuous 
analogue of N 0 /N t . 

To calculate N 0 and P 0 , we use the IP3R model of [7,8], with 
minor modifications described later. 

• /diff =^dift(c/) — c). 



• /serca = V s c"' / (K"' + c"') where K s and n s are obtained from 
[36]. 



• /leak= k\ eak (c s -c). 



• / n includes a basal leak (/leakin), receptor-operated calcium 
channel (ROCC, / r0 cc)j store-operated calcium channel 
(SOCC, / socc ). By using the IP3 concentration (p) as a 
surrogate indicator of MCh concentration, we assume that 
/-occ = V Iotx p. SOCC is modeled by 

/ S occ=K socc < OC c/«occ+^) [13]- 



• J pm =V p c"P/(K;"+c"P). 



Calcium concentration at open channel pore (c p ) does not 
explicidy appear in the equations but is used in the IP3R model 
introduced later. c p is assumed to be proportional to SR Ca 2+ 
concentration (c s ) and is therefore simply modeled by 
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Table 2. Parameter values of the IP 3 R model. 





Parameter 


Value/Units 


Parameter 


Value/Units 


912 


1240 s _1 


921 


88 s - ' 


m 


3s-' 


932 


69 s-' 


926 


10500 s-' 


962 


4010 s-' 


945 


11 s-' 


954 


3330 s-' 


H 


20 s-' 


L 


0.5 s"' 




100 s-' 




100 s-' 




40 s" 1 







doi:1 0.1 371 /journal.pcbi.1 003783.t002 



Cp = Cpo(Cj/100) where c /; o is the value corresponding to 
c s = 100 /(M. Alternatively, c p can also be assumed to be a large 
constant (say greater than 1 00 /iM) without fundamentally altering 
the model dynamics. The choice of c p o is not critical as long as it is 
sufficiently large to play a role in inactivating the open channels. 
All the parameter values are given in Table 1 . 

The data-driven IP 3 R model 

The IP3R model used in our ASMC calcium model is an 
improved version of the Siekmann IP3R model which is a 6-state 
Markov model derived by fitting to the stationary single channel 
data using Markov chain Monte Carlo (MCMC) [5,7,8]. Fig. 1 has 
shown the structure of the IP3R model which is comprised of two 
modes; the drive mode, containing three closed states C\, C2, C3 
and one open state Oe, and the park mode, containing one closed 
state C4 and one open state O5 . The transition rates in each mode 
are constants (shown in Table 2), but 942 and (724 which connect the 
two modes are Ca 2 + -/IP3-dependent and are formulated as 

q 2 4 = a 24 + K24O - W24A24), (4) 



<?42 = «42 + Vnmnh^ (5) 

where m.24, A24, and A42 are Ca 2+ -/IP3-modulated gating 
variables. «24, 042, Via and V42 are either functions of/? or constants 
and are given later. We assume the gating variables obey the 
following differential equation, 



^ = l G (Goi, - G), (G = m 2 4, h 2 4, »?42, ^42), (6) 



where G m is the equilibrium and Xq is the rate at which the 
equilibrium is approached. Those equilibria are functions of Ca 2+ 
concentration at the cytoplasmic side of the IP3R, denoted by c in 
the equations, equal to either c p or Cj depending on the state of the 
channel). They are assumed to be 




Figure 9. Stationary data and fits of q 14 and q 41 . Stationary transition rates of 424 and i]42, <724oo and q42 x , as functions of Ca 2+ concentration 
were estimated and fitted for two [IP3], 1 /iM (A) and 10 /<M (B). Circles and squares represent the means of 924 and 1/42 distributions computed by 
MCMC simulation [7]. Note that MCMC failed to determine the values of q 2 4 and q 42 at c= 1, 10 /iM for 10 /(M IP 3 , as the IP 3 R was almost in the 
drive mode for these cases. The corresponding fitting curves (solid for ^42; dashed for 924) are produced using Eqs. 7-12. 
doi:1 0.1 371 /journal.pcbi.1 003783.g009 
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□ data(p=luM) 
■ data(p=10uM) 




10 10 

[Ca 2 1 fljM) 

Figure 10. Open probability curves for various [IP3]. P„ is equal 
to the sum of probabilities of the IP3R in O5 and Of,. Three 
representative curves correspond to 0.1 /iM, 1 (OA and 10 /(M [IP3] 
(from bottom to top) respectively. Data (average open probability) are 
from [5]. 

doi:10.1371/journal.pcbi.1003783.g010 



'"42 



In 



k 2 

C 2 +k 2 _. 



K -42 

P + k\ 



(8) 



(9) 



(10) 



Hence, we have stationary expressions of q$2 and q24, 

?24oo=«24+' / 24(l-»l24oo^24oo), (H) 



given in Table 2, whereas li H2 is not clearly revealed by 
experimental data. However we have shown that it should be 
relatively large for high c but relatively small for low c for 
reproducing experimental puff data [8]. By introducing two Ca 2+ 
concentrations, C/, and c p> i/, i2 and the state of the IP3R channel 
become highly correlated, so that we can assume 2/, 42 is a relatively 
large value H if the channel is open and is a relatively small value L 
if the channel is closed. Hence, /U 42 is modeled by the logic function 



Mi 



'42 ■ 



H, if the channel is open ; 
L, if the channel is closed . 



Values of L and H are chosen so that simulated Ca 2+ 
oscillations in ASMC are comparable to experimental observa- 
tions. 

The IP3R model reduction 

Here we reduce the 6-state model to a 2-state open/closed 
model. The reduction takes the following steps: 

• The sum of the probabilities of C\, C3 and O5 is less than 0.03 
for any c, so they are either rarely visited by the IP3R or have 
a very short dwell time. This implies they have very little 
contribution to the Ca 2+ dynamics. Therefore, we completely 
remove the three states from the full model. 

• Transition rates of ^26 and qa are about 2 orders larger than 
that of (?24 and q<\2, which allows us to omit the fast transitions 
by taking a quasi-steady state approximation. This change will 
affect two aspects. First, we have Of, = C2q2(,/q€i which allows 
us to combine C2 and Of, to be a new state D, which satisfies 
D= Oeiqa +?26)/?26- Although this means D is a partially 
open state with an open probability of 926 / (?62 + ?26), it can be 
used as an fully open state in the stochastic simulations by 
multiplying the maximum IP3R flux conductance £jpr by a 
factor of (?26/('?62 + ?26)- Secondly, 924 needs to be rescaled by 
qa/iqa + <?26), i.e., the effective closing rate is tf24<762/(<762 + 
qif,)- 

• Due to the combination of C2 and Of,, /./, is accordingly 
modified to 



#42 00=6(42+ V n m n 00 ^42 00 • (12) 

The expressions of as, Ks, ns and ks are chosen as follows so 
that Eq. 1 1 and Eq. 1 2 capture the correct trends of experimental 
values of q24 and 042 (see Fig. 9) and generate relatively smooth 
open probability curves (see Fig. 10), 

V 24 = 62 + 880/(p 2 + 4) fl 2 4 = 1 + 5 l(p 2 + 0.25) 

K 42 = U0p 2 /(j) 2 + 0.0l) a 42 = l-&p 2 /(]> 2 + 0.34) 

^24 = 0.35 £ 42 =0.49 + 0.543p 3 /0' 3 + 64) 

k^24 = 80 £_42 = 0.41 + 25p 3 /(P 3 + 274.6) 



Note that the above formulas are different from the relatively 
complicated formulas used in [8]. The rates, X mit , /l/, 24 and ^ m42 , are 
constants estimated by using dynamic single channel data [4] and 



{H, if the channel is in C4 
L, if the channel is in D (the drive mode) 

Hence, the reduced two-state model contains one "open" state 
D and one closed state C4 with the opening transition rate of 942 
and the closing transition rate of </24(?62/(?62 + <?26)- 

Deterministic formulation of the stochastic model 

Based on the stochastic calcium model and the reduced 2-state 
IP3R model, we construct a deterministic model. We need to 
modify three things that are used in the stochastic model but 
inapplicable to fast simulations of the deterministic model. The 
first is the discrete number of open channels; the second is state- 
dependent use of Cb and c p in calculating 942 and ^24; the last is the 
logic expression of A# 42 . Details of the modifications are as follows, 

• The fraction of open channels (N 0 /Nt) is replaced by open 
probability P 0 which is 70% of the probability of state D. 
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In the stochastic simulations, qiA which only controls the IP3R 
closing is primarily governed by c p , whereas ^42 which controls 
IP3R opening is mainly governed by C/,. Therefore, in the 
deterministic model, we separate the functions of c p and cj by 
assuming m2a, co and /i24oo are functions of c p only whereas 
m&tz, an d IU2a3 ar e functions of Cb only. That is, 
m 2 4oo = «?24a>(t>), h 2 A 02 =h 2 n m (c p ), m 4 2oD ='«42oo(ca) and 
luicc =^42co( c i>)- Here c p = c^ofo/lOO) as defined before. 
To describe an average rate that infinitely many receptors are 
rapidly inhibited by high Ca 2+ concentration but slowly 
restored from Ca 2+ -inhibition. /./, 42 is proposed to be 



full model to a minimal model that still captures the crucial 
features of the full model. First of all, /.„, 42 , X mii and Xh n are 
sufficiently large so that we can assume they instantaneously follow 
their equilibrium functions. Therefore, by taking quasi-steady state 
approximation to /«24, ^24 an d tn^l, we remove the three time- 
dependent variables from the full model. 

By now, the full model has been reduced to a 5D model, 



dc 
dt 



J in J d: 



(21) 



h.,=(\-D)L + DH. 



Based on the above changes, the full deterministic model 
containing 8 ODEs is presented as follows, 



dc 
dt 



- /diff + /leak /serca+/in ^pmi 



deb 
dt 



z 'i>l(JlPR —Jdift), 



dc i=J, -J 
dt 1,1 pm ' 



— = 942(1 -D)-( ; )A 

dt q 62 + 926 



(13) 



(14) 



(15) 



(16) 



dc h 
dt 



= }'l(/lPR — /diff), 



dc, 

~dt 



~T7 =942(1 -D)-( — )D, 



dll42 

dt 



962+926 



, , , , h 42 ). 



(22) 
(23) 

(24) 
(25) 



Second, the rate of change of D approaching its equilibrium, 
= (942962 + 942926 + 924962)/(962 + 925) (calculated from Eq. 24), 
is at least one order larger than those of c, c, and A42, indicating that 
taking the quasi-steady state approximation to Eq. 24 could not 
significantly affect the evolutions of c, c t and A42. That is, 



D- 



942(962+926) 



942962 + 942926 + 924962 



(26) 



dm 2 4 , , c p . 



(17) 



(18) 



We emphasize here that the theory of the quasi-steady state 
approximation has not yet been well established, particularly 
about the rigorous conditions under which such a reduction is 
valid. Thus, our criterion of judging the validity of the reduction is 
checking whether the solutions of the reduced model are capable 
of qualitatively reproducing that of its original model. For this 
model, we find the reduction works. Hence, the full model is 
eventually reduced to a 4D model summarized as follows, 



dlfl42 



C b K 42 



(19) 



dc 
It 



2 J diff ~l~ J. leak J serca ~\~ J'm Jp. 



(27) 



^T=<MtT3 *">' C 20 ) ^=yi(/rPR-/diff), (28) 

at c /j + k_ 42 



where ^24 an d 942 are functions of the gating variables given by 

Eqs. 4 and 5. All the fluxes are the same as those of the stochastic dc, 
model except JwK = kwR(Dq 26 /(q 6 2+q26))(c s -c b ). All the dt 
parameter values of the deterministic model are the same as those 
of the stochastic model and are therefore given in Tables 1 and 2. 



/pm, (29) 



""42 _ < , "--42 / \ /,r,\ 

Reduction of the full deterministic model ~dT~ 42 4+k 3 _ 42 [ ' 

The full deterministic model contains 8 variables which make 
the model difficult to implement and analyze. Thus, we reduce the where D is given by Eq. 26. 
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Inclusion of calcium buffers 

To check the effect of calcium buffers on oscillation frequency, 
we introduce a stationary buffer (no buffer diffusion), as mobile 
buffers are too complicated to be included in the current 
deterministic model. Since we have two different cytoplasmic 
Ca 2+ concentrations, c and c/,, two pools of buffer with the same 
kinetics should be considered. Hence, the inclusion of a stationary 
calcium buffer is modeled by the following system, 



dc 

-j - = -/cliff + Jieak- 



/serca + Jin ~ Jpm ~ k+ (B, - b\)c + k- b\ , (31) 



dt 



■- Yi <7n>R - Jdnt) -k+{B t -b 2 )cb + k_b 2 , 



dc, 
dt 



'- J'm Jm 



(32) 



(33) 



Statistical analysis 

Data analysis is performed on the Ca 2+ traces with relatively 
stable baselines and less noise. A moving average of every 3 data 
points is used to improve the data by smoothing out short-term 
fluctuations (Fig. 2A is an improved result). Due to large variations 
in baseline, amplitude, and level of noise in data, we used two 
thresholds to get samples: a low threshold, 20% of the amplitude of 
the largest spike above the baseline, to initially filter baseline noise 
out; and a relatively high threshold, 50% of the amplitude of the 
largest spike above the baseline, to further remove small spikes that 
cannot initiate waves. For simulated stochastic traces of variable c, 
we first convert it to fluorescence ratio {F/Fq) by using 
F/Fq = c(co + Kd) / (coc + coKd) where the dissociation constant 
of Oregon Green Kd = 0.17 fiM and resting [Ca 2 + ]j c 0 = 0.1 fiM. 
We then used the same sampling procedure mentioned above to 
obtain samples. After samples are chosen, ISIs and spike durations 
are calculated based on the low threshold. Simulated traces used to 
calculate average frequency are about 200-400 seconds long. All 
the samplings and linear least-squares fittings are implemented 
using MATLAB (see Text S3-S4 for Madab codes). 



db\ 
dt 

db 2 
dt 



--k + (B t — bi)c—k-b\, 



= k + (B,-b 2 )c h -k-b 2 , 



(34) 



(35) 



,2+. 



where b (b\ and b 2 ) and B t represent the concentrations of Ca 
bound buffer and total buffer respectively. k + and k- are the rates 
of Ca 2 + -binding and Ca 2 + -dissociation, indicating how fast the 
time scale of the buffer dynamics is. Fast buffer refers to the buffer 
with relatively large k + . In the simulations, we use a fast buffer 
with k + = 100 jiiM^'-s -1 and k- = 100 s -1 and vary B, to test if 
the stochastic model and the deterministic model have a 
qualitatively similar 5,-dependency. Results are given in Fig. 6E. 

Numerical methods and tools for deterministic and 
stochastic simulations 

For the stochastic model, Eqs. 1-3 and ODEs of the four gating 
variables in the IP3R model are solved by the fourth-order Runge- 
Kutta method (RK4) and the stochastic states of IP3R determined 
by the IP3R model are solved by using a hybrid Gillespie method 
with adaptive timing [37]. The maximum time step size is set to be 
either 10~ 4 s (for the 6-state IP 3 R model) or 10~ 3 s (for the 
reduced 2-state IP3R model). All the computations are done with 
MATLAB (The MathWorks, Natick, MA) and the codes are 
provided in Supporting information (Text S1-S2). For the 
deterministic model, we use ode 15s, an ODE solver in MATLAB. 
Accuracy is controlled by setting an absolute tolerance of 10~ 8 
applied to all the variables. 



Supporting Information 

Dataset SI ASMC calcium fluorescence trace data. The 

data files are in Excel format and compressed in a zip file. Each 
Excel file has a name showing their information. For example, 
"S2_SMC6_MCh200nM" means data are from ASMC No. 6 in 
lung slice No. 2 by using 200 nM MCh. In each file, there are four 
columns which represent (from left to right) time(s), fluorescence 
intensity, F/Fo and average F/Frj. 
(ZIP) 

Text SI Matlab code for simulation using 6 state IP3R 
model. 

(DOCX) 

Text S2 Matlab code for simulation using 2 state IP3R 
model. 

(DOCX) 

Text S3 Matlab code for experimental data analysis. 

(DOCX) 

Text S4 Matlab code for simulation analysis. 

(DOCX) 
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